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In this paper we analyze the performance of the Quantum Adiabatic Evolution (QAE) 
algorithm on a variant of Satisfiability problem for an ensemble of random graphs 
parametrized by the ratio of clauses to variables, 7 = M/N . We introduce a set of macro- 
scopic parameters (landscapes) and put forward an ansatz of universality for random bit 
flips. We then formulate the problem of finding the smallest eigenvalue and the excitation 
gap as a statistical mechanics problem. We use the so-called annealing approximation with 
a refinement that a finite set of macroscopic variables ( verses only energy) is used, and are 
able to show the existence of a dynamic threshold 7 = 7 y, beyond which QAE should take 
an exponentially long time to find a solution. We compare the results for extended and 
simplified sets of landscapes and provide numerical evidence in support of our universality 
ansatz. 

PACS numbers: 03.67.Lx,89.70.+c 


I. INTRODUCTION 

An important open question in the field of quantum computing is whether it is possible to 
develop quantum algorithms capable of efficiently solving combinatorial optimization problems 
(COP). In the simplest case the task in a COP is to minimize the energy function E a with the 
domain given by the set of all possible assignments of N binary variables, cr = {oq, .... m v }, 
(jj — ±1. In quantum computation this cost function corresponds to a Hamiltonian H P 
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where the summation is over the 2 N states jcr) forming the computational basis of a quantum 
computer with N qubits. State \o j)j of the j-th qubit is an eigenstate of the Pauli matrix d~ with 
eigenvalue &j. It is clear from the above that the ground state of Hp encodes the solution to the 
COP with cost function E a . In what follows we shall use two equivalent notations for binary 
variables: Ising spins cr. , = ±1 as well as bits Zj = (1 — aj )/ 2 = 0. 1. 

Recently Farhi and coworkers proposed a new family of quantum algorithms for combinato- 
rial optimization that is based on the properties of quantum adiabatic evolution [1,2]. Numerical 
simulations were performed for the study of its performance for satisfiability problems [8]. Im- 
plementation of these algorithms on a quantum computing device is feasible for COPs where the 
energy function £„ possesses a locality property, in a sense that it is given by the sum of terms 
each involving only a relatively small number of bits, that does not scale with N [1, 3, 4], An ex- 
ample of a problem that can have this property is Satisfiability that deals with N binary variables, 
submitted to M constraints, assuming that each constraint involves 0( 1) bits. The task is to find a 
bit assignment that satisfies all the constraints. 

Satisfiability is a basic problem in the so-called NP-complete class [5]. This class contains 
hundreds of the most co mm on computationally hard problems encountered in practice, such as 
constraint satisfaction and graph coloring. NP-complete problems are characterized in the worst 
case by exponential scaling of the run time or memory requirement with the problem size N. A 
special property of the class is that any NP-complete problem can be converted into any other 
NP-complete problem in polynomial time on a classical computer. Therefore, it is sufficient to 
find a deterministic algorithm that can be guaranteed to solve all instances of just one of the NP- 
complete problems within a polynomial time bound. It is widely believed, however, that such an 
algorithm does not exist on a classical computer. Whether it exists on a quantum computer is one 
of the central open questions. 

Running of the quantum adiabatic evolution algorithms (QAA) for several NP-complete prob- 
lems has been simulated on a classical computer using a large number of randomly generated 
problem instances that are believed to be computationally hard for classical algorithms [2, 6, 8]. 
Results of these numerical simulations for relatively small size of the problem instances ( N < 
25) suggest a quadratic scaling law of the ran time of the QAA with N. 

A particularly simple version of Satisfiability is the NP-complete Exact Cover problem that 
was used in [2] to study the performance of QAA. In this problem each constraint is a clause that 
involves a subset of K = 3 binary variables. A given constraint is satisfied if exactly one of its bits 
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equals 1 and the rest of the bits equal 0. In the optimization version of this problem one minimizes 
the energy function E a that is equal to the number of constraints violated by a given bit-assignment 
cr. A generalization of this problem to an arbitrary number K can be called positive I-in-K SAT 
[9]. 

In practice algorithms for NP-complete problems are characterized by a wide range of running 
times, from linear to exponential, depending on the choice of certain control parameters of the 
problem (e.g., in Satisfiability it is the ratio of the number of constraints to the number of vari- 
ables, M/N). Therefore, a practically important alternative to the worst case complexity analysis 
is study of a typical-case behavior of optimization algorithms on ensembles of randomly generated 
problem instances chosen from a given probability distribution. For example, in the case of pos- 
itive 1-in-K SAT one can define a uniform ensemble of random problem instances. An instance 
J consists of M statistically independent clauses, each corresponding to a A-tuple of distinct 
bit-indices uniformally sampled from the interval (1,1V) with probability 1 / (£) . 

In the case of an exponential scaling low for the algorithm’s running times t a it is convenient to 
analyze the distribution of a normalized logarithmic quantity log t a /N. This distribution becomes 
increasingly narrow in the limit of large N where the mean value (log t a ) /N well characterizes the 
typical case exponential complexity of an algorithm. For Satisfiability problem the dependence of 
the asymptotic quantity 

■q = lim (logf 0 )/rV (2) 

h—OG 

on the clause-to- variable ratio 7 = M /N has the qualitative form shown in Fig.l. At some critical 
value 7 = 7 d algorithmic complexity undergoes the dynamical transition from polynomial to 
exponential scaling law. This transition has been studied recently for the case of a variant of 
the classical random-walk algorithm for the Satisfiability problem [10]. Function 77 ( 7 ) is non- 
monotonic in 7 and reaches its maximum at a certain point y c > -/ d . It was discovered some time 
ago [11, 13, 14] that j c is a critical value for the so called satisfiability phase transition: if 7 < y c , 
a randomly drawn instance is satisfiable with high probability, i.e., there exists at least one bit 
assignment cr that satisfies all the constraints (Eg. = 0). For 7 > 7 C instances are almost never 
satisfiable. In the asymptotic limit N — * 00 the proportion of satisfiable instances drops from 1 to 
0 infinitely steeply at 7 = % as shown in Fig. 1. 

The value of ~/ d (unlike y c ) depends on both the problem at hand and the optimization algorithm. 
Comparison of the dynamical thresholds 7 d for different algorithms provides an important relative 
measure of their typical-case performance in a given problem. In this paper we will provide the 
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FIG. 1: Solid line shows the qualitative plot of the normalized quantity r)/r) msx vs M/N (rj m&x is a maxi- 
mum value of rj). Dashed line shows the proportion of satisfiable instances vs M/N. 

analysis of the dynamical threshold for the quantum adiabatic evolution algorithm and also for 
simulated annealing for several versions of the random satisfiability problem. 

II. QUANTUM ADIABATIC EVOLUTION ALGORITHM 

Consider the time-dependent Hamiltonian H{t ) = 7 i{t/T) 

H(t) = (1-t)H b + tH p , (3) 

where r = t/T € (0. 1) is dimensionless “time”, Tip is the “problem” Hamiltonian (1) and Tip 
is a “driver” Hamiltonian, that is designed to cause transitions between the eigenstates of Tip. 
Using dimensionless time and setting ft = 1 the quantum state evolution obeys the equation, 
iTd |^'(r))/(5r = 7f(r)|^(r)). At the initial moment the quantum state l’I'(O)) is prepared to be 
the ground state of 7ft (0) = 7ft g. In the simplest case 

7 i B = -jhai, ^(0))-2-^ 2 ^|cr) : (4) 

j— 1 cr 

where a J x is a Pauli matrix for j-th qubit. Consider the instantaneous eigenstates of 7 i(r) with 
eigenvalues Xk(r) arranged in nondecreasing order at any value of r E (0, 1) 

^(r)|^(r)) = A fe (r)|^(r)) ; (5) 

here k = 0. 1.2, . . . . 2 iV - 1. Provided the value of T (the runtime of the algorithm) is large 
enough and there is a finite gap for all r E (0.1) between the ground and excited state energies, 




Ai(r) — A 0 (t) > 0, the quantum evolution is adiabatic and the state of the system |'£'(t)} stays 
close to an instantaneous ground state, \4>q(t)) (up t0 a phase factor). The state |d 0 (l)) coincides 
with the ground state of the problem Hamiltonian Tip and, therefore, a measurement performed 
on the quantum computer at the final moment t = T (r = 1) will yield one of the solutions of 
COP with large probability. 

The standard criterion for adiabatic evolution is usually formulated in terms of minimum exci- 
tation gap between the ground and first exited states [12] 


T» 




AA ffiin = max [Ai(r) - A 0 (r)] 


( 6 ) 


Here the quantity £ is less than the largest eigenvalue of the operator Tip — Tip [18] and scales 
polynomially with N in the problems we consider. 


III. QUASICLASSICAL APPROXIMATION AND COMBINATORIAL LANDSCAPES 


In the computational basis (1) we have 

H = t E (T \ cr)(a\ - (1 - r) ^<5 [d(cr, cr'), 1] |cr)(</|, (7) 

cr cr ,<y ! 

here 6{m, n) denotes the Kronecker delta-symbol and the summation is over the pairs of spin 
configurations cr and cr' that differ by the orientation of a single spin, d(cr, cr')=l, where 

1 iV 

d(o-,<r r ) =- -o-j-l, (8) 

denotes a so-called Hamming distance between the spin configurations a and cr\ that is the num- 
ber of spins with opposite orientations. Eq. (5) in the computational basis takes form 

A(r)o 0 .(r) = tEo-oJt) - (1 - t) ^ S [ d(<r , cr'), 1] ^(r) (9) 

cr' 

(here we drop the subscript indicating the number of a quantum state in A and In what follows 
we assume that typical energies E a ~ 0(N), but the change in the energy after a single spin flip is 
0(1). This assumption about the energy landscape holds for instances of the Satisfiability problem 
with the clause-to- variable ratio M/N = 0(1), the case of most interest for us (see the discussion 
in Sec. I). 

We now consider a set of functions {Xi = Ciicr.T), l — 1, . . . , /C}, referred to as (com- 
binatorial) landscapes, that depend on a problem instance 1 and project a spin configuration a 
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onto a vector {Xi} with integer- valued components. Prior to considering a specific COP here we 
make certain assumptions about the properties of landscapes and apply them to the analysis of the 
minimum gap in the QAA. 

In particular, we assume that, similar to energy, landscapes {Xi = Ci(cr.X)} are macroscopic 
functions, so that the typical values of Xi are O(N), and possess a certain universality property 
in the asymptotic limit N — > oc. Specifically, the joint distribution of {Ci(cr,T)} over the spin 
configurations a forming the 1 -spin- flip neighborhood of an “ancestor” configuration cr f depends 
on a problem instance I and spin configuration a f only via the set of parameters {X' t — Ci(a^X)}. 
We then define a quantity 

1 ^ 

| {*;}) = - ao) 

d(<r.cr / )=l k~l 

x; = a (a-'.i), 

In effect, the above universality property of landscapes implies that the set of all possible 
spin configurations a is divided into “boxes” with coordinates {Aj} where A; — Ci{er), and 
P ({X;} | {X/}) (10) represents the transition probability from box { Xi } to box {A/}. In particu- 
lar, it obeys Bayes’ rale 

p({x t } | {xi}) n({x{}) = p({xd | {a,}) Q{{x ( }), (ii) 

where Q({A;}) is the number of different spin configurations in the box {A)}. 

We consider energy to be a smooth function of landscapes 

E^ — E ({A;}) , Xt^Qia,!), (12) 

so that \dE/dXi\ — 0(1). Furthermore, we assume that, on one hand, the change in C'i(cr.T) 
after flipping one spin is 0(1), for typical problem instances. On the other hand, we assume 
that correlation properties in a neighborhood of a box {A;} described by P ({A;} | {A'/}) vary 
smoothly with box coordinates on a scale 1 < \SX (( <C N. Therefore if we write the transition 
probability in the form 

P ({A/} | {A)}) = p ({XI - AJ: {x,» , = Xi/N}, (13) 

then p({ki }: {xi}) is a steep function of its first argument: it decays rapidly in the range 1 < 
\ki\ <C N for each (-component. However this is a smooth function of its second argument: it 
varies slightly when coordinates xi change on a scale |c)a; ; | -C 1. 
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One can show that under the above assumptions the quantum amplitudes 6 a corresponding to 
the smallest eigenvalue depend on the spin configuration a only via the coordinates of thes box 
{A'j} to which it belongs. Then we look for the solution of (9) in the following form: 






{X, = G(<t, I)}. 


(14) 


where \(p({Xi} , t)\ 2 gives the probability of finding the system in the box {A';}. Plugging (14) 
into (9) and making use of (1 1),(12) we obtain: 


X(t)(p(X, t) = tE(X)(p(X. t) — (1 — t)N ^ L(X. X'}) <p(X', r), • (15) 

X' 

X = {X 1 ,X 2 ,...,X IC }, (16) 


(hereafter we use the above shorthand notation for the set of landscapes). In (15) we introduced 

L 


L(X,X') = L(X',X) = P(X'j X) (17) 

P(X) = 2~ N Q(X). 


where P(X ) is a probability that a randomly sampled configuration a belongs to a box X. We 
shall look for a solution of (15) in the WKB-like form 

ip(X,T) - exp (—W(X, r)) , (18) 


so that 

A (r) = tE(X) - (1 - T)N^2L(X,X')e W{X - r) - w{X '’ T) . 

X' 

We now introduce scaled variables (c f. (13)) 



_y_ 

7n’ 


(19) 


( 20 ) 


and also 

w(x,r) = lw(X,T), c(x) = jfE(X), a(x) = T logn(X), (21) 

where s(x) is an entropy function. Based on (17) and the properties of the transition probability 
(see Eq. (13) and discussion after it) we assume that the sum over X' in (19) is dominated by terms 
with |X' - X| = 0{l). Then we can use an approximation 


W (X'. r) - W(X, t) k Vw • (X' - X) + C(l/N), 


( 22 ) 
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where Vw = dw(x, Fj/cbc. Plugging (22) into (19) and making use of Eqs. (13), (17), (20) and 
(21) we obtain after some transformations: 

g = h(x, Vw; T), (23) 

h(x. p; T) = g(x) - F y^p(k;x)e~ k (v ^ 2 ~ i ' p) . 

k 

(here Vs = ds(x)/dx). This is a Hamilton-Jacobi equation for an auxiliary mechanical system 
with coordinates x, momenta p = Vw, action w, Hamiltonian function h(x, p: F) and energy g. 
Using the symmetry relation 

p(k; x)e~ k VS//2 = p(—k; x)e k Vs ^ 2 , (24) 

that follows directly from Eqs. (11) and (17) we obtain that the minimum of w(x, F) over x where 
Vw = 0 necessary corresponds to the minimum of the functional: 

/(x,T) =e(x) - H?(x), (25) 

where f(x, F) = h(x , 0, F) and 

£(x) =p (Vs/2: x), p(y: x) = p (k; x) e _k y . (26) 

k 

The summation in (23) and (26) is over components k[ of k in the range ki € (— oo, oo). In what 
follows, we shall refer to p( y; x) in (26) as a “Laplace transform” of p(k; x). 

We note that £(x) = ^ x/ L(X / .X) and one can use Bayes rule and inequality of Cauchy- 
Bunyakovsky in (17) to show that that the positive-valued function £(x) is bounded from above, 

0 < f’(x) < L This shows that the analysis of the effective potential based on the WKB approxi- 
mation (22) is self-consistent in the asymptotic limit N — ► oo. 

It follows from the above analysis that the ground-state wavefunction y>(x, T) = y?(X,r) 
is concentrated in x-space near the bottom of the “effective potential” given by the functional 
/(x,T), i.e. near the point x ¥ (T) where /(x.F) reaches its minimum. In this region S « |x r Ax, 
where matrix A is positive definate, and according to (18), the wavefunction has a Gaussian form 
with the width oc i/Vn. 

The ground-state energy g = g(T) is given by the value of the effective potential / (25) at its 
minimum 

.?(r) = /(x„(r),r) t (2?) 

a/(x, F)/5x| x=jt _ (r; = 0, /(x, F) > g(T). 
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We note that as F — 0 the shape of the effective potential /(x, F) approaches that of the energy 
function s(x) and therefore its minimum x„(F) — *• x 0 where x 0 is a minimum of c(x). It can be 
shown that in this limit the ground-state eigenvalue approaches the minimum energy value s(x 0 ) 
and the eigenvalues of A -1 approach zero (and so does the characteristic width of the wavepackage 
c(x, F)). The spin configurations that belong to a box x 0 in x-space correspond to the solutions of 
the optimization problem at hand. It is clear that one of the solutions can be recovered with high 
probability after a measurement is performed at the end of the “quantum annealing” procedure. 

Variational Ansatz: For cases in which the set of macroscopic variables (A';} is not suf- 

ficient (in statistical sense (13)) to describe the dynamics of the quantum algorithm, one can 
still implement the above procedure as an approximation, using a variational method. Intro- 
ducing a Lagrangian multiplier A, one looks for the minimum of the functional F(<p, A) = 
(®\H\d>) — \((<p\<p) — 1), using a variational ansatz (14) for the wavefunction. The solution of 
the variational problem is provided by Eqs. (18)-(27). The smallest eigenvalue g (27) corresponds 
to the value of the Lagrange multiplier at the extremum, A = tN g, and the maximum of the 
variational wavefunction corresponds to the minimum of the effective potential / (25). 

A. Global bifurcations of the effective potential 

However, in the case of a global bifurcation where the effective potential /(x, F) possesses 
degenerate or nearly degenerate global minima, the answer is modified. If for some value of 
F = F,, a global bifurcation occurs, in our example this would mean that for this value of F, two 
values of x, x+ and x~ give a global m inim um to /(x, F). In such a case, the smallest eigenvalue 
is not doubly degenerate; rather an exponentially small gap AA min between the ground and first 
excited state is developed, itself being proportional to the overlap between two wave-functions, 
peaked around x+ and xF respectively. 

To estimate the overlap we note that at T, the two global minima of the effective potential 
/(x, F„) correspond to the two coexisting fixed points of the Hamiltonian function in (23) with 
zero momentum and the same values of energy g, 

df/dx = dh/dx - dh/dp — 0 
x = x~, p = pf=0 ; 5 (x,p;F.) = gt = g~ . 


( 28 ) 

( 29 ) 



10 


Then to logarithmic accuracy we have 

~ log A^min = J dt' [x(f')P(0 - h(x(t ') . p(t’)) } + 0(1/N), (30) 

where ( x(t), p(£) ) is a heteroclinic trajectory connecting the two fixed points of (23) 

x(t) = dh/dp, p (t) = - dh/dx , (31) 

x(f — > ±oo) = x*, p(f — »• ±oo) - 0 . 

From the algorithmic perspective this means that when F gets close to F,, it has to change 
exponentially slowly (cf. Sec. II and Eq. ( 6 )). This could be called a critical slowing down in 
the vicinity of a quantum phase transition. If simulated annealing (SA) is used and a similar 
phenomenon occurs, the value of the temperature T» is the point where a global bifurcation occurs 
in the free energy functional 

/(x,T) = s(x) -Ts(x). (32) 

By comparing the free energy functional (32) with the functional (25) corresponding to “quantum 
annealing” (QA), we note that in QA the quantities F and £(x) play the roles of temperature and 
entropy in (SA), respectively. 

We note in passing that a similar picture for the onset of global bifurcation that can lead to 
the failure of QA and (or) SA was proposed in [18, 19] for the case where the energy E a is a 
non-monotonic function of a single landscape parameter, a total spin °r 1 ° ^is case the 
dynamics of QA can be described in terms of one-dimensional effective potential [20, 23]. 

IV. THE MODELS 

An instance of a Satisfiability problem with N binary variables committed to M = 7 N con- 
straints (where each constraint is a clause involving K variables) can be defined by the specifi- 
cation of the following two objects. One of them is an M x N matrix Q, the rows of the matrix 
are independent A’-tuples of distinct bit indexes sampled from the interval ( 1 . N). The m-th row 
of Q defines the subset of the K binary variables involved in the m-th clause. The second object 
is a set of boolean functions B = {b m }, with each function encoding a corresponding constraint. 
A function b m = b m [a Smi , a g ^ , , <J SmK ] is defined over the set of 2 K possible assignments 
of the string of K binary variables involved in the m-th clause. The function returns value 1 for 
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assignments of binary variables that satisfy the constraint and 0 for bit assignments that violate it. 
Then the energy function equals to the number of violated constraints 

.v/ 

E<J — E <?@) = M - T i a e ml ■ a Cm2 ’•••’ a s mK j • 

m= 0 

here X = (Q, B) denotes an instance of a problem. 

The matrix Q defines a hypergraph Q that is made up of the set of N vertices (corresponding 
to the variables in the problem) and a set of M hyperedges (corresponding to the constraints of 
the problem), each one connecting K vertices. An ensemble of disorder configurations of the 
hypergraph corresponds to all the possible ways one can place M = jN hyperedges among N 
vertices where each hyperedges carries K vertices. Under the uniformity ansatz all configurations 
of disorder are sampled with equal probabilities (i.e., rows of the matrix Q are independently and 
uniformly sampled in the ( 1 ,N) interval). 

Boolean functions b m may also be generated at random for each constraint with an example 
being random K-SAT problem [16, 17]. However here we consider slightly different versions of 
the random Satisfiability problem that are still defined on a random hypergraph Q but have a non- 
random boolean function b m = b, identical for all the clauses in a problem. One of the problems 
is Positive 1-in-K Sat in which a constrain is satisfied if and only if exactly one bit is equal 1 and 
the other K- 1 bits are equal 0. The boolean function b for this problem takes the form 

' K 1 — a 

b{ai . oti = 6 E 1 — , 1 (Positive 1-in-K Sat). (34) 

.p= i 

a p = ±1, p = 1,2. .... K. 

We shall also consider another problem. Positive K-NAE-Sat. in which a clause is satisfied unless 
all variables that appear in a clause are equal (”K-Not-All-Equal-Sat”). The boolean function b 
for this problem takes the form 

K . , 

b[ai, q 2 , . . . , Oat] = 1 — ^ S ^ T — , 0 (Positive K-NAE-Sat). (35) 

S=± 1 _p=l 

Both problems are NP-complete (Appendix A). It will be shown below that they are characterized 
by the same set of landscape functions. 
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VI LANDSCAPES : ANNEALING APPROXIMATION 


For a particular spin (cr) and disorder (Q)) configurations, all clauses can be divided into 2 A 
distinct groups according to the values of the binary variables that appear in a clause. We will label 
the different types of clauses by vectorial index a = {ai, ... ,ax}, a p = ±1. We now divide 
the set of 2 N spin configurations into boxes identified by certain numbers of clauses of each type, 
N M a , and also by the Ising spin in a configuration N q 


1 M K 

M a = M a (cr, Q) = - IT 5 [° Smp > a p] > 

m= 1 p— 1 

l jV 

q = q(tr) = 


(36) 

(37) 


Different boxes correspond to macroscopic states defined by the set of parameters (q, {M a }) with 
q £ (—1, 1) and J2 a M a — 7- The energy function can be expressed via (36) as follows (cf. 
(33H35)): 


K 

€({**.})= T-'E&Mn, M. 

m=0 




K 

K — 2m. '^2 a p f 

p = 1 


(38) 


where the form of the coefficients depends on the problem: 


Cm 


5{m, 1] (Positive 1-in-K Sat) 

1 — S[m , 0] — 5[m, K] (Positive K-NAE-Sat) 


(39) 


In the following we compute an approximation to the effective potential (25), using the land- 
scape functions (36), (37). According to (26) it depends on the entropy function s(q, {M a }) and 
the transition probability (13) between different macroscopic states. Recalling that variables q and 
M a are normalized by the factor N we study the probability of transition, p(n. {r Q }; q, {M a }), 
from the state (q, {M a }) to the state (q + n/N, {M a + r a /N}). The Laplace transform of p with 
respect to n, {r Q } has the form (cf. (26)) 


p(0,{yay,q : {M a })= J2 e- en - E “ Vara P(n,{rahq,{M a }), (40) 

n, (r a } 

We assume that all binary variables are also subdivided into distinct groups based on their value 
cr — dil and a vector k with integer coefficients k p a indicating the number of times a variable 
appears in a clause of type a in position p. Clearly, consistency requires that k p a = 0 unless 
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a p — a. We now define a quantity c a ^ which is equal to the fraction of spins with given cr, k. For a 
spin configuration cr there exists a set of coefficients {c^k} with elements of the set corresponding 
to all possible values of a and k (there will be many 0’s in a set for each spin configuration). In 
general, there axe exponentially many sets {c^k} that correspond to a macroscopic state (g, {M a }) 


Y^° c *M = <h y>£c, ik = iV/ Q (p = 0, 1, . . . , K). (41) 

o\k tr,k 

Coefficients {c a k} are concentrations of spin variables with different types of “neighborhoods”. 
We shall assume that in the limit of large N the distribution of coefficients c G y corresponding to 
the same macroscopic state (41) is sharply peaked around their mean values (with the width of the 
distribution cc A r “ lj/2 ). 

Under the above assumption we can immediately compute the Laplace-transformed transition 
probability (40) in terms of the coefficients c Gi k- Indeed, consider flipping a spin with value cr 
and neighborhood type given by vector k. This will change the total spin by — 2a and for each 
clause of type a and index p E (l,i\) the value of NM a will decrease by k On the other 
hand, for the clause type ol = a(p y a) obtained by flipping a bit in p-th position in a, NM& is 
correspondingly increased by k Hence the Laplace-transformed transition probability is 


P (0 : q, {-Wq}) = ex P 29cr + Y2( y °~ y & (p’<*)) k a 

a, k L 

where the coefficients c^k are set to their mean values in a macroscopic state (41)). 


(42) 


A. Entropy and coefficients c a k in a macroscopic state defined by q and {A/ a } 

Here we use the annealing approximation to estimate the mean values of c^k and also of a 
macroscopic state (g, M a ). We start by introducing the concept of annealed entropy. Let j\f be 
the number of spin configurations subject to some constraints. In general, it is a function of the 
disorder realization. The annealed entropy is defined as the logarithm of its disorder average: 
Sann = In (A/)* Note that for the correct, quenched, entropy the order of taking a logarithm and 
disorder average is reversed. 

Since in the random hypergraph model all disorder configurations are equally probable, an- 
nealed entropy is given as s ann = lnyV^g — InA/g, where Ms,g is the total number of spin and 
disorder configurations and Mg is the number of disorder configurations. 
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For enumerating all possible disorder configurations we depart slightly from the traditional 
random hypergraph model. In our model all clauses are ordered (two disorder configurations where 
any two clauses are permuted are deemed different), clauses can be repeated (the same clause can 
appear twice), the order of variables in a clause is important (two disorder configurations are 
different if the order of variables in any clause is changed), and finally, variables can be repeated 
in a single clause. This change does not alter the underlying physics, since the probability that 
two identical clauses appear is infinitesimal, and a variable enters a clause twice in at most 0 ( 1 ) 
clauses, which can be safely neglected. As regards the distinction between the disorders with 
permuted clauses, this only introduces a combinatorial factor which cancels out. The advantage is 
that each disorder can be represented as a sequence of M A^-tuples of integers from 1 to A r . 

We will first compute the annealed entropy of a macroscopic state (q } {M a }) under addi- 
tional constraints: we fix the values c ffj k and compute the annealed entropy as a function of 
q. {A / a } 3 {c ffi k }. Recalling that M a are the numbers of clauses of a given type scaled by N, 
and the total number of clauses is 7 JV, we obtain the number of joint spin-disorder configurations 
as a product of the following factors: 

(i) the number of ways to assign types to clauses (Nj)\/ ^(AMq)-^ 

(ii) the number of ways to assign types to variables N\/ 


(iii) for all p,a, the number of ways to permute the appearance of variables in p-th position of 
clauses of type a: (NM a )]/Yl a ^( 


Consequently, the annealed entropy is given by 


-f- ( K ~ 1) ^2 A/ a I* 1 A dot + 7 In 7 — 7 A'. 

a 

(43) 

In the large N limit we replace c^k by their annealed averages, i.e., the values that maximize the 
annealed entropy. In its simplest form, we place no constraints on c a% k except consistency require- 
ments (41). Associating Lagrange multipliers A and In n p a with these constraints, the expression 
for the entropy can be rewritten as 


•®ann [{*-cr,k}! Q- { ^ ^ Or, ,k ltl 

C7,k 


C<rM 


n^ ! ) 


^ann [$: {A/qc}J 


min ( -A q + Y] M a In ^ + In Z[ A, {n p a } j 
Y, M a In M a + 7 In 7 - 7 K. 


( 44 ) 
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The values of e CT .k are given by 


c,,=^n^f~/k> a L 


p, Of 


and Z is given by 

Z = exp (x + z 1 ]^') + ex P (~ x + Z Z 6 '- a P- ~ 


a p 


a p 


The values of the Lagrange multipliers A, p? a are related to q, {M a } via 

31nZ 


dX 

dhiZ 


= q, 


^ dfj, p a Ma ' 


From here we obtain the expression for the Lagrange multiplier /i£ 

Met _ 1 t OLp q 


Then introducing a new notation 




\r^ 1 + a p P 

= 2 _ ~ 2 ~^ 
p, Of 

M± = 14 

V-> 01 


(45) 

(46) 

(47) 

(48) 

(49) 


we obtain 


Z = eV+ +e"V-. 


Md 


2M± 

T±~q 


(50) 


(51) 


Then the entropy can be rewritten in the following form 


5ann[<?, {-V7q}] = ~ Xq -j- M+ In -i- M_ In Tin M a In M a 4 - 7 In 7 - 7 AT (52) 


1-q 


1 + 9 - V- =z 1 ~ q 


We now use the following equations 

eV + = Z . - _ _ 0 

2 ' 2 

and obtain the expression for the second Lagrange multiplier A 

1+9, 1 + 9 1-9, 1-9 

■ In 


-Xq = -- 


2 ln L_!_i nZ + jK. 


(53) 


(54) 


2 2 

Upon substitution of A from the above into the expression for s ann (52) we finally obtain the 
annealed entropy 

•Sam, [ 9 > {-W Q }] - - q tanlU 1 q - In — ^ — + M+ In + M - In 

- M a In M a + 7 In 7 . (55) 

Of 

Also the coefficients c CTj k are given by (45),(46) with Lagrange multipliers given in (49) and (54). 
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B. Effective potential 


Consider a factor £(x) = (df/dF) s (25), (26) in the expression (25) for effective potential with 
x = (q, {M a }). It follows from (26) that to find this factor we need to evaluate the Laplace- 
transformed probability (40,42)) at 

9 — ^dSznn/dq, Hot — ^ds^aa/dMoc- (56) 

This is where the Lagrange multipliers come in handy as we can immediately claim that 


8q A 


dSa.nn/dMa = E ln “F _ ln M, 


(57) 

(58) 


(J-c 


Note that in differentiating with respect to above we omitted the constant term. This is permis- 
sible since only differences dq^/dM* - dq^/dM^ appear in Eq. (42). A further refinement is 
to write 




1 ' gpg = K In g ■ + S2 cr p tanh 1 q. 

9 2 4— ' v 


Using this in the Eqs. (26), (42), we obtain 


£(q, {M a }) = ± E II ^ Ep/(V 

tr,k p,ot \ 


) tanh 1 q 


Mr 


Jt£ 


/*&!• 


Since \ )T .(cy - cr' p ,) = a p (where a' is obtained from a by flipping, p-th bit) and also 


M q // q = tanh -1 1 


the expression is considerably simplified 


£(q, {Af a }) = Jexp 


= E 


where the sum is over pairs (a, a') that differ in exactly one position 




K 


o E _ a 'j \ = L 


To evaluate Z we write 


\A - 7 


p - 1 


Ve^-e^- 


2 ( M+ 

exp 


V 1 - ? 


M- 


1 +q 1-7 


(59) 


(60) 


( 61 ) 
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and the expression for £ becomes 


£{q- Woe}) = y/l - g 2 exp 


2 E< Q .a'> v%-V Q - A/ + A/_ 


(62) 


Vi -<? 2 1 + ? 

here A/_ are given in (50). 

We note that the effective potential f(q. {M a }) = c({A/ Q }) — Tf(g. {M a }) is symmetric with 
respect to permutation of individual components in {A/ Q } corresponding to different orders of 
-l’s and +l’s in the vectorial index a. We look for the minimum of f(q. {M a }) using symmetric 
ansatz 


m 


-l 


K 


M a =W 1 M m , rn = Y J ] —^ 2 

p=l 


(63) 


where m is the number of -l’s in a. Substituting (63) into (62) and rewriteing 

" 2 Em=o \A m ± 1 )(- K ~ ~ ™)M m M m +i 


Z( q: {M m }) = VTT 


q- exp 


W ± j Em=o ( K ~ 2m)Ad r 

l -? 2 




(64) 


where we defined £ [q, {M m }) = {M a }). The effective potential is then 

f(q- {Mm}) = s {{M m }) - T£ (q, {M m }) ( QA ), 


(65) 


with energy given in (38). In the case of the SA algorithm the corresponding free-energy functional 
(32) is 

/(<?■ {Mm}) = z{{M m }) ~ Ts(q . {A4 m }) (5.4). ( 66 ) 

where the entropy function equals 

s(q.{M m }) = - q tanli -1 q + (yK - 1) In V l — - 

( K \ K V4 

- i Y^( K ~ 2m )M m J tanh^ q -^2 MmlnW^-. (67) 

\m=0 / 771—0 \m) 

If we were to use an even smaller set of macroscopic parameters (e.g. only the energy e) 
we can still employ formula (64) with the proviso that unspecified variables should be taken to 
equal their most likely values, i.e. those that maximize the entropy s(q, {M m }) not the landscape 
£(q t {M m }). For example, in the case of energy-only landscapes, £ = £(s), the values q. {M m } 
that max imiz e s{q.{M m }) for a given energy s and number of hyperedges jN (X(m=o M m = 7 ) 
should be computed and then substituted into the expression for £ (64). 
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1-in-K 

K-NAE 

K 

Id 

7c 

Id 

7c 

n 

J 

- 

0.805 

- 

2.41 

4 

0.676 

0.676 

- 

5.19 

5 

0.557 

0.609 

- 

10.7 

6 

0.475 

0.548 

19.8 

21.8 

7 

0.416 

0.500 

34.9 

44.0 

8 

0.371 

0.461 

61.7 

88.4 

9 

0.335 

0.428 

109 

177 

10 

0.305 

0.400 

196 

355 


TABLE I: Annealing bounds for dynamic (7^) and static (7 C ) transition for positive 1-in-K SAT and positive 
K-NAE SAT for different values of the number of variables in a clause K. 


We compute, within the annealing approximation, the point of static transition (cf. Fig.l), 
where the entropy of the macroscopic state with zero energy vanishes, 5(0) = 0, and the dynamic 
transition 7^; for connectivities 7 > 7 d an effective potential (65) exhibits a global bifurcation for 
some T = T*. The resulting values are given in Table I ( see also Figs. 2 and 3). 



FIG. 2: Static j c (circles) and dynamic 7 d 

(crosses) transition for positive 1-in-K SAT vs 
K. 


FIG. 3: Static 7 C (circles) and dynamic j d 

(pluses) transition for positive K-NAE SAT vs 
K. 


In Fig. 4 we plot time variations of the landscape parameters, M m = A4* m , corresponding 
to the global minimum of the effective potential. In Fig. 5 we plot a time-variation of the scaled 
ground-state energy g given by the value of the effective potential at its minimum. Singular be- 
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havior corresponding to the first-order quantum phase transition at certain r = t„ (F = F») can 
be clearly seen from the figures. Plots in Figs. 4 and 5 correspond to precisely the static transition 
7 — 7 C for the case of K = 4 in 1-in-K SAT problem. In the region jd < 7 < there are 



FIG. 4: Plots of the landscape parameters M m = M* m at the global minimum of the effective potential, 
vs r for K = 4 (1-in-K SAT problem). Curves labelled 0-4 correspond to 0/7 through M, 4 / 7 . 



FIG. 5 : Scaled energy of adiabatic ground state g vs t 
for K=4 (1-in-K SAT problem). 

exponential (in N) number of solutions to Satisfiability problem but the runtime of the quantum 
adiabatic algorithm to find any of them also scales exponentially with N. This is a hard region for 
this algorithm. We note, that in the limit of K — > 00 the annealing approximation becomes exact. 
Together with the fact that for large K '/ d and y c seem to be distinctly different provides evidence 
that this result (existence of hard region for quantum adiabatic algorithm) is robust. 
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FIG. 6: Results of numerical simulations and their comparison with theory. Depicted are Laplace transforms 
of Mi for l-in-3 SAT. Numerical results: curves that have different colors correspond to different random 
problem instances; curves of same color correspond to different random bit strings. The dashed black line is 
a theoretical result based on the annealing approximation. The insets (a)-(e) depict instances with 10°, 10 4 , 
10°, 10 6 , and 10' binary variables. Since the error is not recognizable we replot in (f) a magnified section 
of inset (e). The bit strings were sampled with q = 0.422, Mo = 0.048, M\ = 0.416, M 2 = 0.123, 
= 0.013, corresponding to M/N = 0.6. These values correspond to the energy E^/2 and they are 
shifted by 10% from the most likely values of q, {M m } for this energy (this shift is » N l > 12 ). We also note 
that for l-in-3 SAT numerical simulations give static phase transition at 7c « 0.62). 








21 


VI. UNIVERSALITY PROPERTY FOR TRANSITION PROBABILITIES 


Here we study the universal features of the transition probability in (10) for the set of macro- 
scopic variables corresponding to the (normalized) total Ising spin q and numbers of clauses of 
different types {M m } (38) (the type of a clause is equal to the number of unit bits involved in the 
clause). For simplicity, we shall focus in this section on the case K = 3 only. 

To clarify the above choice of macroscopic variables we consider an auxiliary quantity: a 
conditional probability distribution of the macroscopic variables (q, {M m }) over the set of all 
possible configurations cr obtained by flipping r bits of the configuration cr'. The first moments of 
this distribution corresponding to A4 m , 

Mm = O) y) s [d{cr' - cr), r] M m (cr, J), m = 0,...,K, (68) 


can be easily computed by counting the number of ways one can flip r bits in configuration cr' to 
transform a K - bit clause of m' type (i.e., with m! unit bits) into a clause of the m-th type 


p.m'—O 


p j \m — nr 


A r - K 

r — 2p — m -r ml 


(here we use the convention Q) = 0 for m < 0 and nn > n). In the double sum above values of 
A4 m ' are multiplied by the number of possible ways to flip three groups of bits: p unit bits in a 
clause of ra'-type, p + m~m f zero bits of this clause, and r — 2p — m + m! bits of the configuration 
c 7 f that do not belong to the clause. Similarly, one can show that the first moment corresponding to 
the variable q equals q'( 1 — 2 r/N). It is clear that dependence of the first moments on the ancestor 
configuration cr' is only via the variables q\ M' m for that configuration. 

In the limit, r > 1, the above conditional distribution has a Gaussian form with respect to 
q and M m . Elements of the covariance matrix {cr 9 ) = C?(r), and correspondingly, the 
characteristic width of the distribution is 0(r 1//2 ). For a configuration cr' randomly sampled in the 
box (q, {M m }) the r.m.s. deviation of the elements of (cr') from their mean values in the box 
is 0(N l,/2 ). It is clear that in the limit r > A rl/2 the covariance matrix elements can be replaced 
by their mean values for the macroscopic state (<?, { Ai m })• Therefore in this limit the conditional 
distribution after r spin flips starting from some macroscopic state depends only on the values of 
(q : in this state (universality property). 

One can show that for r < N l/2 the conditional distribution after r spin flips can be expressed 
via the distribution (10) with r — 1, using a standard convolution rule. However for r = 1 the 



form of the distribution is non Gaussian and we were not able to establish universality properties 
in the general form. Instead we performed a series of numerical studies. In Figure 6 we present 
the results of numerical simulations and the comparison with analytic results within the annealing 
approximation. One can see that the theory is in very good agreement with experiment. 


VII. CONCLUSION 

We have formulated an ansatz of landscapes and studied the complexity of the quantum adia- 
batic algorithm within the annealing approximation and found the existence of a dynamic transi- 
tion and a hard(exponential) region above that dynamic transition. However, a similar analysis of 
simulated annealing did not reveal any phase transitions. We explain this as follows. The anneal- 
ing approximation should fail for sufficiently small energies. It is commonly known that simulated 
annealing can find suboptimal solutions with very small energies very efficiently, but it takes an 
exponentially long time to actually reach the ground state. The annealing approximation does not 
correctly describe very small energies and cannot be used to establish its complexity. Note that 
we can reconcile this with the fact that the annealing approximation becomes exact in the limit 
K — * co: if the annealing approximation fails for E < Ek we expect that E K is decreasing to 
zero as K increases. However for any finite K , the free energy computed within the annealing ap- 
proximation is free from any singularities indicative of a phase transition. To study the complexity 
of simulated annealing one needs to use the tools of spin glass theory, in particular, the replica 
trick [15-17]. 

In contrast, in our analysis of the quantum adiabatic algorithm, we observed a first-order phase 
transition, and, importantly, it happens for energies E ~ 0{E oo ) (where E ^ is the expected 
energy at infinite temperature E^ — E z . Moreover, the energies on both sides of the 

transition, relative to E ^ seem not to change appreciably with increasing K. Since the annealing 
approximation for this range of energies can be used, the prediction for the dynamic transition 
should survive, though the exact numerical values may acquire corrections. We have recomputed 
the dynamic transition with simplified energy-only landscapes (see Fig. 7). For 1-in-K SAT one 
can clearly see that the relative correction quickly diminishes. We believe that same happens for K- 
NAE SAT if sufficiently large K ’ s are considered. If this indeed holds, it serves as a corroboration 
that our results are correct numerically for large K. The idea of using energy-only landscapes was 
present in [7] as well as [21] and [22]. A jump in the time-dependence of the expected energy 
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value was seen in numerical simulations [8], indicative of first-order phase transition, though a 
different ensemble was considered (only instances having a unique solution were considered). 



4 5 6 7 B 9 10 


K 

FIG. 7 : Relative difference between predictions for the dynamical phase transition point in the case of full 
( 7 4 ) and energy-only (-yf) landscapes vs of K for 1-in-K SAT (crosses) and K-NAE SAT (pluses). 

We emphasize that the annealing approximation employed in this paper essentially neglects 
fluctuations due to disorder, and describes the transition as a global bifurcation between two 
macroscopic states (pure states) and the complexity is due to tunnelling between them. In con- 
trast, spin glass theory predicts the existence of an infinite number of pure states[15]. Secondly, 
affirming our results for large K ignores the structure of the problem, since that limit corresponds 
to the so-called random energy model, where one does not expect to do better then 0(2 A/2 ) via 
any quantum algorithm. Consequently, the complexity could be determined not by the unique 
minimum gap, but by a cascade of level repulsion. Numerical studies, however, support the pic- 
ture with a unique minimum gap. Also, the first-order phase transition occurs for large energies. 
Although it is absent for small K, we believe that a better approach (as compared to annealing 
approximation) will reveal it. Moreover, we believe that the order of the transition will remain 
unchanged, suggesting that the disorder may be irrelevant for the determination of the order of the 
phase transition and, consequently, for the complexity of the quantum adiabatic algorithm. That 
is, the exponential complexity is not due to the true combinatorial complexity of the underlying 
random optimization problem but rather due to peculiarities of the driver term and a particular 
ensemble of random instances considered. In fact, for a symmetrized variant of the exact cover 
problem, the same phenomenon was observed - the exponential slowdown - although the problem 
did not possess any randomness [18. 19]. In fact, a ground state of that problem could be found 
in O(N) time. However, it was possible to modify the driver term in the annealing Hamiltonian 
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[20, 23] to circumvent the slowdown. It is quite possible that a similar change of driver term can 
achieve same goals in present case, although we have not analyzed this scenario. In such a case, 
one would have to go beyond the annealing approximation to study the complexity. 
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APPENDIX A: ON THE NP-COMPLETENESS OF POSITIVE 1-IN- A" SAT AND 

POSITIVE A'-NAE-SAT. 

We set out to prove that both Positive l-in-A' Sat and A'-NAE-Sat are NP-complete. It is 
straightforward to see that it takes a polynomial time to verify the assignment, so these problems 
are in NP. We now prove that they are as hard as the Satisfiability problem, which is NP-complete, 
by showing that any boolean formula can be represented as an instance of these. 


1. Positive l-in-A' Sat 

A clause of type (x. . . . , x, y) necessarily implies x = 0 and y = 1; hence we can represent 
constants 0 and 1. A clause of type (0, .... 0. x. y) implies x = ~y. Finally, a clause of type 
(0, . . . , 0. x, y, z) is equivalent to a 3-clause (x, y, z) so that we can restrict ourselves to A" = 3 
without losing generality. 

For K = 3, immediately observe that three clauses (x. z, u')(y, z, u")(u, u' . u") with free vari- 
ables u, u /, u” implies z = -<(x A y). This basic building block is in fact sufficient to build any 
boolean formula, as a result, any boolean formula can be cast as an 1-in-AT SAT formula. 
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2. Positive A' -NA E-Sat 

A clause of type (x, . . . , x, y) necessarily implies x = -i y, and (x, . . . , x, ?/. z) is equivalent 
to (x.y.z) so we once again restrict ourselves to K — 3. In contrast to 1-in- A'' problem, we 
shall require a non-trivial representation of false or true . We will use pairs of variables to denote 
variables of the boolean formula. Pairs 00 or 11 will represent value false and pairs 01 or 10 will 
represent true. 

The next building block, (x, y, t)(y. z, t)(z, x y t) ensures that t = 1 if the majority of x, y, z 
are 0 and t == 0 if the majority are 1. We shall use a shorthand f(t: x,y 5 z) to denote this. The 
expression f(zi\ x l3 y 1: y 2 )f(z 2 : x 2 , j/i, 2/2) then ensures z = x A y where x, y, z are represented 
as pairs xix 2) yiV 2 , z x z 2 as indicated above. The operation of negation is trivial to represent: if 
x = X 1 X 2 then ->x = (“ixi)x 2 . These two are sufficient to construct any boolean formula. 


APPENDIX B: NEXT ORDER APPROXIMATION FOR LANDSCAPES 


A better approximation for the values of critical clause-to- variable rations can be obtained if 
we specify the constraint that the distribution of vertex degrees be Poisson (as it is supposed to be 
in a random hypergraph [24]). To be precise, we specify that 


/W C(j '. k n ^ ( zL/ ^ c {k P ) — 

<r,k p \ oc / 


iva 


Consequently, with this constraint the following expression for is obtained: 


C cr , k ~ C {kp}~ 


z {h) 


where we use 


z {k P } = e> n f ~ ) + e a n f E) 5 ( a ? + ^ 


(Bl) 


(B2) 


(B3) 


Annealed entropy can be rewritten in the form 


■Sann[<7' {M a}] - min < -A q 4- V' M a In - + In Z[ A, i - A/ a In 4- 7 In 7 - 7K, 

l Pltt ^ J 


where In Z is given by 


InZ = V^c {/Cp} ln Z {M . 

{kp} 
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The equations relating q. {M a } and A. {ji p a } are given in (47),(48). 

Similarly to Sec. VB we will use the notation (50). Since In Z depends only on A and / a± , 
d In Z/dp? a depends only on <j p . Therefore, M a / n p a = A/ CTp //v„- Correspondingly, 

s ann [q. {M Q }1 = min ( — A q + M+ In — + M_ In — + In Z } - Y M a In M a 4- 7 In 7 - 
A, Mi [ M+ M- J „ 

For convenience, we introduce new variables 


Mi 


= /j, e 


±h 


(B4) 


We then readily obtain 


In Z = 7 In n + Y Ck In [2 cosh(A + kh)}, 

k 

(where k = k p and c k — and drops out of the expression for s ann altogether: 

s ann [q, {A/ a }] = min | —A q — (M + — M-)h + Y^ c k In [2 cosh(A + kh)] | (B5) 

- Y M <* ln Ma + M + 111 M + + M ~ ln M ~ + 7 ln 7 “ 1 K - 


It is easy to see from this expression what the equations for A, h are: 

'Y c k tanh[A + kh] — q, 

k 

y kc k tanh[A + kh] = M + — A/_ . 


(B6) 


We now turn our attention to the function £(q, {A/ Q }) given by (42) with 6 and y a evaluated from 
Eqs. (56).(58). The computation of e Va ~ yn ’' yields 

a p 

1 u , 

(B7) 


/AL. 


e y«-y a > = t ftE) 

V-V-VMi/ V A/a 


Multiplied by pP a this becomes 


ha e Va y “ = nrr " M- 


y/MTMl 

The expression for P,(q, {A/ a }) can be written in the form (cf. (59)) 


(B8) 


<(s, {«,}) = e z n 

{fc p } oM 




p L 
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with the internal sum running over k consistent with a set of {A; p } (Bl). Substituting the quantities 


defined above this becomes 


{*M) ~~ V y C {kp} 

{k p } 

After some transformations we finally obtain 


ZM P /*£«*(<*-*) 




-<a.a'> V ° 


cosh [A + A;/i] 


(BIO) 


where A, h are given by (50),(B6). Using symmetric ansatz (63) it is straightforward to calculate 
from (BIO) the restricted function £(q, {M m }) (cf. (64)). We must note however, that although 
this represents a next-order improvement over annealing approximation, the relative changes in 
and computed with this improved approximation are nearly imperceptible (~ 10 -4 ). 


